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ABSTRACT 

Context. In protoplanetary disks, dust grains coagulate with each other and grow to form aggregates. As these aggregates grow by 
coagulation, their filling factor (p decreases down to <c 1 . However, comets, the remnants of these early planetesimals, have ~ 0. 1 . 
Thus, static compression of porous dust aggregates is important in planetesimal formation. However, the static compression strength 
has been investigated only for relatively high density aggregates (0 > 0.1). 
Aims. We investigate and find the compression strength of highly porous aggregates {(p <g: 1). 

Methods. We perfom three dimensional A'-body simulations of aggregate compression with a particle-particle interaction model. We 
introduce a new method of static compression: the periodic boundary condition is adopted and the boundaries move with low speed 
to get closer. The dust aggregate is compressed uniformly and isotropically by themselves over the periodic boundaries. 
Results. We empirically derive a formula of the compression strength of highly porous aggregates (0 « 1). We check the validity of 
the compression strength formula for wide ranges of numerical parameters, such as the size of initial aggregates, the boundaiy speed, 
the normal damping force, and material. We also compare our results to the previous studies of static compression in the relatively 
high density region (<p > 0.1) and confirm that our results consistently connect to those in the high density region. The compression 
strength formula is also derived analytically. 

Key words, planets and satellites: formation - methods: numerical and analytical - protoplanetary disks 



^ 1. Introduction 



Planetesimal for mation is a key issue in planet formation in pro- 
toplanetary disks (iHav ashi et al.ll 19851 IWeidenschilling & Cuzzil 
[l993i) . However, collisional growth of dust from sub-micron 
sized dust to kilo-meter sized planetesimals is still unknown. 

In the growth process, one of the most important but unre- 
solved problems is the internal structure evolution of dust aggre- 
gates. Dust internal structure is important in planetesimal forma- 
tion because dynamics of dust aggregates in protoplanetary disks 
is determined by coupling between gas and dust, in other words, 
size and internal density of dust aggregates. In the early stage of 
dust coagulation in protoplanetary disks, the collision ene rgy of 
the aggregates is too low to cause collisional compression ([Blum 



interactions (iDominik & Tielenslll997l: IWada et al.l[2007i |2008L 
' ' " " kl2008l[2009[ 



2009l:ISuyama et alj2008ll2012tlPaszun & Dominiki2008il200S 
Seizinger et al .I l2012h . 



npr 

2004 l iOrmel et al.ll2007HZsom et al.ll2010il20T ItlOkuzumi et al 



20121) . As a result, the internal mass density p decreases down 
to p < 1.0 g cm"^. 

Both theoretical and experimental studies have shown that 
mutual collisions lead dust aggregates to have their frac- 
tal dimension D ~ 2, whic h is so-called b allistic c luster- 
cluster a ggregation (BCCA) (ISmimoyI Il990t iMeakinI 
Kempf e talll 19991: iBlum & Wurml l2000l: iKrause & BlumI 



1991 



20041 : 



Paszun & Domini fl l2006h . The dust aggregates would be grad 



ually compacted or disrupted in coagulation because of the in- 
crease in impact energy. Such compaction has been investigated 
with numerical A^-body simulations considering particle-particle 



In most of previous studies investigating dust growth in pro- 
toplanetary disks, dust grains have been assumed to have con- 
stant internal mass d ensity for sirnplicity (iNak agawa et al II 198 lb 

lUni 



iTanaka et al.1 120051: iBrauer et al.l 120081: Mmstiel et al.r i2010l) . 
However, dust porosity evolves during dust growth in protoplan- 
etary disks in reality. In recent dust coagulation calculations, 
porosity evolution has been considered based on experimental 
and theoretical results ( Ormel et ani2007l: lOkuzumi et al 1 120091 
I2OI2I : IZsom et alJl201 ll) . Thev also suggested that p decreases 
as p <g: O.lg cm"^. 

In the most recent work, though dust grains have size distri- 
bution, the dominant coagulation mode has shown to be s imilar- 
size collisions of dust aggregates (lOkuzumi et al ] 120121) . As a 
result, their fractal dimension is approximately equal to 2 and 
their internal mass density p has shown to become 10"^ g cm"-' 
(equivalent to be the filling factor (p - 10"^ for ice particles with 
a density of 1 .0 g cm"-'). Such flulfy dust aggregates are believed 
to become planetesimals. Since comets in our solar system, 
which would be remnants of planetesimals, have their internal 
mass density of ~ 0.1 g cm"-' (A'Hearn 2011), dust aggregates 
must be compressed from p «: 0.1 g cm"-' top ~ 0.1 g cm"-' in 
protoplanetary disks. 
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Compression at dust aggregate collisions has been inves- 
tigated in previous studies. When colUsional impact energy 
exceeds the critical ener gy , dust aggregates are compacted 
by th eir collsion (e.g. iD ominik & TielensI [T997t ISuvama et al.l 
120081: IWada et al.1 [200^12008. ,20091). However, the colhsional 
compressi on is not effective to compress dust aggregates plan- 
etesimals dOkuzumi et al ] l2012h . 

One of the other compression mechanisms in protoplanetary 
disks is static compression by disk gas or self-gravity. The static 
compression strength of dust aggregates has been investigated 
both experirnentally and numerical ly (Paszun & Dominik 2008; 
iGuttler et al]|2009l: TSeizinger et al.ir2012l) . However, they exam- 
ined only relatively compact aggregates with p > O.lgcm"^ 
because their initial aggregates are ballistic particle-cluster ag- 
gregation (BPCA) clusters. Because p decreases down to p « 
O.lg cm""* at least in the early stage of dust growth, we need to 
reveal the static compression strength with p <*: O.lg cm"-'. 

In this work, we investigate static compression of highly 
porous aggregates with p < O.lg cm"^ by means of numerical 
simulations and analytical approach. It is challenging to per- 
form numerical simulations of static and uniform compression of 
highly porous aggregates. Because such porous aggregates have 
low sound speed, we have to compress them in much slower ve- 
locity than in the case of compact aggregates, as will be shown 
in our simulations. Such a slow compression of the fluffy aggre- 
gates costs much computational time. 

In previous numerical studies of static compression, a dust 
aggregate is compressed by a wall m oving in one direction 
(Paszun & Dominik 2008; S eizinger et al. 2012). However, this 
method has disadvantages to reproduce uniform and isotropic 
compression. There are also side walls which do not move. 
Such side walls also obstruct the tangential motion of monomers 
in contact with the walls, causing artificial stress on the aggre- 
gate and restructures them. Moreover, since they measure the 
pressure with the force on the moving wall, the side walls may 
affect the pressure measurement. In the present work, we de- 
velop a new method to reproduce static compression. Instead of 
the walls, we adopt periodic boundary condition and the bound- 
aries are getting closer to each other With this slowly-moving 
periodic boundaries, the aggregate is compressed uniformly and 
naturally. The periodic boundary condition also enables us to 
represent a much larger aggregate than that inside the computa- 
tional region. This saves the computational time remarkably. 

This paper is organized as follows: we describe the model 
of our numerical simulations in Section |2] We show the results 
of our simulations and find the compression strength in Section 
[3] We confirm the obtained static compression strength formula 
analytically in Section |4] We present our conclusion in Section 



setting in this section. In our simulations, the aggregate is grad- 
ually compressed by its copies over the moving periodic bound- 
aries. This is an appropriate method to simulate uniform and 
isotropic compression. We also describe the boundary condi- 
tion in this section. Since we do not have walls to measure the 
pressure in the periodic boundary condition, we use a similar 
manner of pressure measurement in molecular dynamics simu- 
lations. We also introduce the method of pressure measurement 
below. 

2.1. Interaction Model 

We calculate the direct interaction of each connection of parti- 
cles, taking into account a ll me chanical interactio ns modeled by 
iDominik & TielensI (Il997h and lWada et aP (l2007h . The material 
parameters are the monomerradius ro, surface energy y. Young's 
modulus E, Poisson's ratio v, and the material density po. Table 
[T]lists the values of the material parameters for ice and silicate. 

We perform A^-body simulations with ice particles except 
for one case with silicate particles. In protoplanetary disks, ice 
particles are the most dominant dust material beyond the snow- 
line. Moreover, computational time required for calculation of 
ice particles is less than that of silicate. Thus, we adopt ice par- 
ticles in the most simulations. We also perform a silic ate case to 
compare with a previous study (ISeizinger et al.ll20T2l) . 

The critical displacement still has a discrepancy between the- 
oretical ((cjit - 2 A) and experimental (fp Ht - 32 A) studies 
(IDominik & Tielens"1997';'He im et al.ll999l) . We adopt the same 
parameter of Wada et al.. (.20111) . ^cnt = 8 A as a typical length 
for i ce particles, and ^crit = 20 A for silicate particles to compare 
wifli ISeizinger et al.l d2012h . 

^crit is related to strength of rolUng motion. The rolling mo- 
tion between monomers is crucial in compression. The rolling 
energy Erow is the energy required to rotate a particle around a 
connecting point by 90 °. The rolling energy can be written as 



(1) 



(see lWada et al.l ( |2007|) for details). In the case of ice monomers, 
for example, iSroU = 4.37 x 10"' erg for ^ciit = 8 A. 

We use normalized unit of time in our simulations. In case 
of ice particles, the normalized unit of time is 



/ ^1/2 7/6 N 



to = 0.95 



'0 



1/6 



= 1.93 X 10" 



(2) 



which is a characteristic time and represents approximately the 
oscillation time of particles in contact at the critical collision ve- 
locity ( see lWada et al.1 (|2007|) for details). 



2. Simulation Setting 



2.2. Damping force in normal direction 



We perform three dimensional numerical simulations of com- 
pression of a dust aggregate consisting of a number of spherical 
monomers. As the initial aggregate, we adopt a BCCA cluster 
We solve interactions between all monomers in contact in each 



The normal force between two monomers is repulsive when the 
monomers are close or attractive when they are stretched out. 
Thus, normal oscillations occur at each connection. For reaUstic 
particles, such oscillations would dissip ate due to viscoelastic- 
e p. Interactions between mo nomers in contact are formu- ity or hysteresis in the no rmal force (e.g. iGreenwood & JohnsonI 
lated bv lDominik & TielensI (Il997h and reformulated with using 120061: iTanaka et"ani2012h . For such damping of normal oscilla- 

(120071) . We use the interaction tion, we add an artificial normal damping force to the particle 

interaction model, following the previous studies (Su vama et al.l 
2008; Paszun & Dominik 2008; Seizinger et al. 2012). 

Assuming that two particles in contact have their position 
vectors Xi and X2, respectively, the contact unit vector iic is de- 



potential energies b y IWada et al. 



model proposed by IWada et al.l (120071) in this work. We briefly 
sum marize the particle interaction model and material constants 
(see IWada et al.l (|2007|) for details). Moreover, we describe the 
additional damping force in normal direction and the simulation 
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Table 1. Material parameters in our simulation 



Material 


ice 


silicate 






( Salllc dS oclZlIliicr cl al. KZaJlZ.)) 


Monomer radius ro [jum] 


0.1 


0.6 


Surface energy y [mJ m"^] 


100 


20 


Young's modulus E [GPa] 


7.0 


2.65 


Poisson's ratio v 


0.25 


0.17 


Material density po [g cm"-'] 


1.0 


2.65 


critical rolling displacement ^cnt [A] 


8 


20 



fined as 
rir 



(3) 

\X\ -X2\ 

(see Figure 2 in IWada et al.1 (l2007h '). We introduce a damping 
force between contact particles in normal direction, defined as 



^ damp — '^n 

to 



(4) 



where k„ is the damping coefficient in normal direction and niQ is 
the monomer mass. The adopted value of k„ is an order of 0.01. 
To show that the result is independent of the normal oscillation 
damping, we perform A^-body simulations with the damping fac- 
tor k„ as a parameter. 

The timescale of damping is 



"^damp 



(5) 



for k„ - 0.01, is much shorter than the simulation timescale, 
which is typically ~ lO^fo. We show that the obtained com- 
pression strength is independent of the artificial normal damping 
force in our simulations (see Section lj4t . 



2.3. Uniform Compression by Moving Boundaries 

We adopt the periodic boundary condition in our simulations. 
The aggregate in the computational region is surrounded by its 
copies, as shown in Figure[T] Initially, we set a cubic box whose 
sides are periodic boundaries with a size of L to be larger than 
the aggregate. Thus, the initial BCCA cluster is detached from 
its neighboring copies over the periodic boundaries. In our sim- 
ulations, we gradually move the boundaries to the center of the 
aggregate to get closer to each other. As a result, the aggregate 
sticks to the neighboring copies and is compressed by them in a 
natural way. Therefore, the aggregate in the computational re- 
gion corresponds to a small part of a whole large aggregate. In 
other words, although the number of particles in numerical sim- 
ulations are limited because of computational cost, the periodic 
boundary condition enables us to investigate a large aggregate, 
such as a ~ cm-sized dust aggregate in protoplanetary disks. 

Another advantage of the periodic boundary condition is that 
we do not need to introduce the wall for compression. In the 
previous A^-body simulations of static compression, dust aggre- 
gates ar££om£ressed_byu^ingt^ aggre- 
gate dPaszun & Dominikll200alSeizinger et al.ll20I2h . The wall 
itself may have some artificial effects on such experiments. For 
example, the wall moves in one direction and thus this may be 
different from isotropic compression. Besides, wall-particle in- 
teraction is different from particle-particle interaction, and thus 
it must be treated carefully. In contrast, the periodic boundary 



i» 

































Fig. 1. Schematic drawing of the periodic boundary condition. Each 
of the box illustrates a boundary box with a side length L for all direc- 
tion. When the boundary starts to get closer, the aggregate sticks to the 
neighboring aggregates over the boundary and compressed by them. It 
should be noted that this picture is illustrated in 2D direction, but our 
simulations are performed in 3D. 



condition does not need walls for compression because a dust ag- 
gregate is compressed by the neighboring aggregate over the pe- 
riodic boundary. In addition, the periodic boundaries in three di- 
rections make it possible to compress the aggregate isotropically. 
Note that we calculate not only the interactions of particles in 
contact inside the computational region but also the interactions 
of the particles in contact across the periodic boundaries. Thus, 
no special treatment of interactions, which is wall-particle inter- 
actions in the case of simulations with walls in previous studies, 
is required when a particle crosses the periodic boundaries. 

The computational cubic region has length L and the coor- 
dinates in x,y, and z directions are set to be -L/2 < x < L/2, 
-L/2 < y < L/2, and -L/2 < z < L/2, respectively. We adopt 
periodic boundary conditions for all directions to reproduce a 
part of a large aggregate. L decreases with time t, L - L{f). The 
initial size of the box Lq is adopted as the maximum size of the 
dust aggregate in x, y, and z directions. 

With the settings above, we move the boundaries of the com- 
putational region toward the center of the region. The velocity at 
the boundary is given by 



Vb = -— L(f), 



(6) 
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where Cy is a dimensionless parameter (we call Cy the strain 
rate parameter hereafter). Owing to this definition of the bound- 
ary speed, the aggregate is compressed at a constant strain rate 
independent of the region scale L. 

The box size decreases with the constant rate Cy in three 
directions. This corresponds to isotropic compression. Since 
^ = 2vfe, the box size is written as 



L-Lq exp -2Cv — 



(7) 



Therefore, the whole time of compression is order of to/Cy. Typ- 
ically we chose Cy - 3 x 10"^ and thus the compression time is 
~ 3 X lO^fo ~ 0.6 ms. 

When a particle crosses a periodic boundary, the velocity 
should be treated carefully to reproduce the quasi-static com- 
pression with periodic boundary condition. Figure |2] illustrates 
how to calculate the velocity of particles across the periodic 
boundary. When a particle goes out of the computational re- 

1 time step 











(Vx 




Vb 




Vb 






— > 


^< — 


L 





(vx+2vbyy) 



periodic boundary 



Fig. 2. Schematic drawing to illustrate how the particle velocity is 
calculated when a particle crosses a periodic boundary. For simplic- 
ity, we consider this situation in two dimensional field but we actually 
calculate this in three dimensional situation. We consider that a dust 
particle is close to the boundary in the left figure. In the next time step, 
the particle crosses the boundary (dashed circle in the right figure). We 
put the particle on the other side of the boundary as expressed in Equa- 
tions ([8} and l llOb . The velocity component is converted as expressed 
in Equations ^ and dllb . This treatment well reproduces the isotropic 
compression in the velocity field. 

gion across the boundary at x = L/2, we relocate the particle to 
the opposite side (i.e., from the boundary at x = -L/2). In that 
case, the position of the particle in x direction is converted as 

XI— »x-L (8) 

Since the two boundaries at x = -L/2 and x - L/2 have a rel- 
ative velocity of 2vb, the x-component of the velocity v t of the 
particle is also converted as 



Vx + 2vb. 



(9) 



Owing to the conversion of v^, the velocity of particle against 
the boundary which the particle crosses does not change before 
and after the crossing. For a particle across the boundary at x = 
-L/2, the position and the velocity are converted as 



X I — > X + L 

v., I — > Vx - 2vt. 
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We also have the same treatments for particles across the bound- 
aries aty - +L/2 and z - +L/2. 

We introduce the constant strain rate at the boundaries for 
scaleless discussion. However, the initial aggregate is not mov- 
ing. As the simulation starts, if all the particles in the aggre- 
gate are not moving, only the particles close to the boundaries 
have initial velocity. This is not a constant strain rate. To re- 
produce the scaleless constant strain rate initially, therefore, we 
initially give all monomers the velocity smoothly connected to 
the boundary speed. The initial velocity is expressed as 



v(r) = Vb X 



Lo/2' 



(12) 



where r is the position vector of the monomers. 



2.4. Pressure Measurement 

In previous studies, a dust aggregate is enclosed by walls and 
the pressure was calculated by measuring the force exerted on 
the walls by the dust aggregate. In this work, a dust aggregate 
is compressed by themselves because of the periodic boundary 
condition. Therefore, we introduce another method to measure 
the pressure on the aggregate. We calculate the pressure of 
the dust aggregate with the standard way in molecular dynam- 
ics si mulations using the virial theorem as follows (e.g., Hail^ 
[T99I . 

Let us consider a virtual box that encloses an aggregate under 
consideration. We define the force acting from the walls of the 
virtual box on the particle / as W,, and the sum of the forces from 
other particles on the particle / as f , . The equation of motion of 
the particle ; is given by 



(13) 



We take a scalar product of both sides of the equation with r, and 
take a long time average of the both sides with time interval r. 
The left-hand side becomes 



1 r d^Vi 1 

m- r, ■ -— - = m- 

T Jo dfi T 



dri 
It 



1 r dr, dri , 

-m- I — -dt. (14) 

T Jo dt dt 



The first term in the right-hand side vanishes in the limit of 
T oa. We define taking a long time average in t as (), Taking 
summation of all particles and a long time average of Equation 
( HJI ), we obtain 



(15) 



The first term in the right-hand side is related to pressure P. The 
pressure is an average of all forces acting on the wall from all 
particles. Using the normal vector n of the wall surface directed 
outward, the force received by the wall which has an area dS is 
PndS . Therefore, 



2 



r, ■ Wi 



Pn ■ rdS 



-3PV. 



(16) 



(10) This equation is obtained by taking surface integral as 



(11) X-^'' = //^^^'^ = X(i-^l"l)'^^'^- '''' 
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The translational kinetic energy K, averaged over a long time, is 
given by 



Using K and P, Equation (fTsT l gives an expression of P as 
P^^KIV +UY,ri-F\ IV. (19) 



(18) 



We define the force from particle j on particle / as /,j The force 
Fi can be written as a summation of the force from another par- 
ticle as 



Fi — ^ fij- 



(20) 



J*' 



Using fi j - -fjj, we finally obtain the pressure measuring for- 
mula as 



P^^K/V+Uj](n-rj)-f,}i IV. 



(21) 



The first term in right-hand side of the equation represents the 
translational kinetic energy per unit volume and the second term 
represents the summation of the force acting at all connections 
per unit volume. This expression is useful to measure the pres- 
sure of a dust aggregate under compression. We do not need 
to put any artificial object such as walls in simulations because 
Equation (l2Tl i is totally expressed in terms of the summation of 
the physical quantities of each particle, which are the mass, the 
position, the velocity, and the force acting on the particle. In our 
calculations, we take an average of pressure for every 10,000 
time steps, correspondent to 1000 fo because we set 0. 1 fo as one 
time step in our simulation. 

As mentioned in Section l2!2l the adopted damping force cor- 
responds to rapid damping of normal oscillations. Thus, the 
kinetic energy of random motion rapidly dissipates. This cor- 
responds to the static compression and thus the compression 
strength is determined by the second term of Equation (|211 . 

3. Results 

The top three panels of Figure |3]show snapshots of the evolution 
of an aggregate under compression in the case where = 16384, 
Cv = 3 X 10"^, = 0, and ^crit = 8 A. The top three pan- 
els have the same scale but different time epochs, which are t 
= 0, 1 X lO'^fo, and 2 x lO^'fo, respectively. The white parti- 
cles are inside the computational region enclosed by the periodic 
boundaries, while the yellow particles are in the neighbor copy 
regions. (For visualization, we do not draw particles in the front- 
and back-side copy regions.) The bottom three panels represent 
the projected positions onto two-dimensional plane for the cor- 
respondent top three figures. We confirm that the dust aggregate 
is compressed by their copies from all directions. As the com- 
pression proceeds, the aggregate of white particles is compressed 
by the neighboring aggregate of yellow particles. We focus on 
how high pressure is generated by quasi-static compression in 
numerical simulations. Our numerical simulations have several 
parameters; the size of the initial BCCA cluster, the compres- 
sion rate, the normal damping force, and the critical displace- 
ment (corresponds to the rolling energy). We investigate the de- 
pendence of the pressure on these parameters, by performing a 
lot of runs with different parameter sets. Although we assume 
ice aggregates in most runs, we also investigate cases of silicate 
aggregates to compare them with previous studies. 



3.1. Fiducial run: obtaining ttie compression strength 

We put a BCCA cluster as the initial aggregate. The BCCA clus- 
ter is created by sticking the copy of the aggregate from random 
direction. The results depend on the random number of the ini- 
tial condition, which is the shape of the BCCA aggregate. To 
avoid the dependence, we take arithmetic averages of ten simu- 
lations of different initial conditions. The pressure is measured 
using Equation (|2T] ) at each run. We define the filling factor of 
an aggregate as 



VqN 
V ' 



(22) 



where Vq is the monomer volume, is the number of monomers 
of the aggregate, and V is the volume enclosed by the bound- 
aries, which has a length of L. The filling factor also can be 
written (p - pIpQ. Figure |4] shows that the measured pres- 
sure as a function of the filling factor (p{t). The parameters of 
the simulations are A? = 16384, Cy = 3 x 10"^ /tn = 0.01, 
and ^crit = 8 A. The corresponding E^oW is 4.74 x 10"*^ erg for 
^crit = 8 A. Each colored line in Figure HI a) shows each sim- 
ulation with the different initial shape of the aggregate. Figure 
Htb) shows the arithmetic average of the pressure measured in 
ten different runs. Each line shows in different ranges of 0. The 
lowest is determined with the largest size of the initial bound- 
ary boxes of the ten runs. We find that the compression strength 
is well reproduced by 



(23) 



where Pq - 4.74 x 10^ Pa. We analytically discuss why the com- 
pression strength is proportional to (fr' in Section H] In the high 
density region (0 > 10"'), the measured strength deviates from 
the line of P = Pq(P^ . This is because the dissipation mechanism 
changes in the high density region (see Section U4l i. The devia- 
tion in the low density region (0 3 x 10"^) is partly caused by 
a finite boundary speed (or compression rate) as discussed in the 
next subsection. Another reason of the deviation in the low den- 
sity region is related to the density of the initial BCCA cluster. 
The filling factor of BCCA 0bcca is estimated as. 



»BCCA 



Vbcca 



3/2 



A^" 



-1/2 



(24) 



where we use the radius and the volume of a BCCA clus- 
ter, rsccA ^ V573A^'^Vn and V RrcA = (47r/3)r^pp^, respec- 
tively (e.g.. ISuvama et al.l 120081) . For A^ = 16384, we obtain 
(pBCCA ~ 3 X lO""'. In the early stage of compression, (p is 
lower than (pBCCA because the initial BCCA clusters are apart 
from each other This space between BCCA clusters would also 
cause the deviation from the line of P = Pq(P^ . 

Now, we discuss the coefficient Pq of the compression 
strength. Wada et al. (2008) shows that £,-011 is important in the 
collisional compression strength. Thus, £1011 is expected to be 
also important in the static compression strength. Considering 
that the characteristic volume is monomer's volume ~ rjj, we 

suppose Po = £^roii/'"o, based on dimension analysis. Therefore, 
the compression strength can be written as 



(25) 



We analytically discuss and confirm this equation in Section |4] 
We also plot this equation in Figure Efb). This figure clearly 
shows that the result is well fitted by Equation ( 
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t=0 (ct>=0.0003) t=l X 10*to (ct)=0.002) t=2 X 10% ((1)=0.01) 




-20 -15 -10 -5 S 10 15 20 ' -20 -15 -10 -5 5 10 15 20 ' -20 -15 -10 -5 5 10 15 20 



(Mm) (urn) (urn) 

Fig. 3. Snapshots of the evolution of an aggregate under compression in the case of N = 16384. The top tfiree figures are three dimensional 
visualization. They have the same scale with different time epoch. The white particles are inside a box enclosed by the periodic boundaries. The 
yellow particles are in neighboring boxes to the box of white particles. For visualization, we do not draw the copies in back and front side of the 
boundaries but only 8 copies of the white particles across the boundaries. Each bottom figure represents projected positions onto two-dimensional 
plane of all particles in each corresponding top figure. The gray points in the bottom figures correspond to the positions of the white particles in 
the top figures and the yellow points con'espond to those of the yellow particles in the top figures. Scales are in yum. 



We show that compression strength is proportional to ^crit, 
that is proportional to the rolling energy iSroU in Section [33] We 
also confirm that Equation (l25T l is applicable to the case of dif- 
ferent ro in the silicate case. 



Vd, can be written as 



Vd = |2vb| = 2- 



.Cv 

to 



1/3 



(26) 



3.2. Dependence on the boundary speed 

To statically compress the aggregate, we should move the bound- 
ary at a sufficiently low velocity not to create inhomogeneous 
structure. Figure |5] shows the dependency on the strain rate pa- 
rameter. Each line shows the average of ten runs. The fixed 
parameters are = 16384, k„ = 0.01, and ^crit = 8 A. The 
strain rate parameter Cy is equal to 1 x 10"^, 3 x 10"^, 1 x 10"^, 
3 X 10"^, and 1 x 10"^, respectively. The higher Cy, the higher 
pressure in the low density region is required for compression. 
This is mainly caused by the ram pressure from the boundaries 
with high speed. 

When the compression proceeds and the density becomes 
higher to reach the line of Equation (l25t . the pressure follows 
the equation. From Figure |5] Cy = 3 x 10"^ creates sufficiently 
low boundary speed. The boundary speed can be calculated as 
a function of 0. Using Equation ^ and - {4/3)7TrQN/L^, the 
velocity difference between a boundary and the next boundary. 



In the case of Cy = 3 x 10 ^, Vd = 12.7, 5.9, and 2.7 cm/s for 
- 10"^, 10"^, and 10"', respectively. 

Here, we discuss the velocity difference of boundaries, com- 
paring with the effective sound speed of the aggregates. The 
effective sound speed can be estimated as 



Cs.eif 




£'roll P 



mo 



(27) 



where we use Equation ( l25T l. Using the rolling energy of ice 
particles, Cs,eff is given by 



Cs.eff ~ 1.1 X lO"*!^ cm/s. 



(28) 



Therefore, in the case of Cy = 3 x 10"^, v^i is not sufficiently 
low in the beginning of the simulation, where the aggregate has 
a low filling factor. However, the boundary velocity difference 
reaches lower than the effective sound speed when (p > IQ^~. 
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Fig. 4. (a) Pressure P in [Pa] against filling factor (f>. The ten thin solid lines show the results for the initial BCCA clusters with different initial 
random numbers and thick solid line shows the arithmetic average of the ten runs, (b) Pressure P in [Pa] against filling factor (p. Same as the thick 
solid line in (a) plotted with a dotted line of Equation J25t . The parameters are A' = 16384, Cy = 3 x 10"^, k„ = 0.01, and ^^it = 8 A. 




Fig. 5. Pressure P in [Pa] against filling factor (p with different strain 
rate parameter Cy. Each line shows the average of ten runs of the fixed 
strain rate: Cv = 1 x 10"', 3x 10"^ 1 x 10"*, 3x 10"^ 1 x 10"^ The other 
parameters are the same for every ten runs : A' = 16384, k„ = 0.01, and 
^crit = 8 A. The dashed line is Equation l l25b . 

3.3. Dependence on the size of the initial BCCA cluster 

To confirm that Equation (l25l l is valid in the lower density re- 
gion, we perform the simulations with the different number of 
particles, which is equivalent to the different sizes of the ini- 
tial dust aggregates. Figure |6] shows dependence on the number 
of particles of the initial BCCA cluster. The initial numbers of 
particles are 1024, 4096, and 16384. The other parameters are 



Fig. 6. Pressure P in [Pa] against filling factor with different number 
of particles A^. Each line shows the average of ten runs of the fixed 
number of particles: A' = 1024, 4096, and 16384. The other parameters 
are Cy = 3 x 10"', = 0.01, and ifciit = 8 A in the case of N = 
1024, 4096, and Cy = 1 x 10"', k„ = 0.01, and = 8 A in the case of 
A' = 16384. The dashed line is Equation l l25b . 



Cv = 3 X 10"', /tn = 0.01, and ^crit = 8 A in the case of = 1024 
and N = 4096, and Cy = 1 x 10"', kn = 0.01, and ^.rft = 8 A 
in the case of N = 16384. We chose lower Cy in the case of 
= 16384 in order to investigate the strength in lower (f> region. 
Each line represents the average of ten runs for each simulation 
as in Figures |4jb) and |5] We draw the averaged line from the 
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lower (j) than that in Figure |5] In such a low (p region, we con- 
sider for some runs that the pressure is zero because the aggre- 
gate is isolated from the copies of the aggregate over the periodic 
boundaries. Except for the initial deviation in low (f>, all lines 
have a good agreement with Equation (IZST i where (p ^ 0.1. The 
result has the good agreement in lower (p for runs with larger A^. 
Therefore, we conclude that the formula Equation (l25t is valid 
for0 < 0.1. 

3.4. Dependence on the normal damping force 

As described in Section 12.21 we adopt the norm al damping 
force t o reduce the normal oscillations in addition to IWada et aU 
(l2007h . To confirm that this damping factor does not affect the 
simulation results, we set the damping factor as a parameter 
Figure Q shows dependence of pressure on the normal damping 
factor kn. The fixed parameters are = 16384 Cy = 3 x 10"^, 




Fig. 7. Pressure P in [Pa] against filling factor (fi with different nor- 
mal damping force. We put the same ten initial conditions varying the 
normal damping force with = Q, ka = 10"^, and = 10'. Each 
line shows the result of one run. The other parameters aie N = 16384, 
Cv = 3x 10-\and^,,i, = 8 A. 

and ^crit = 8 A. Each line represents the result of one run for 
ka = 0, 10"^, and lO', respectively. This figure clearly shows 
that the normal damping force does not affect the simulation re- 
sults. 

As mentioned in Section lTTl the compression strength in the 
low density region (0 < 0.1) is expected to be determined by 
the rolling motion. In order to confirm this, we calculate the 
total energy dissipations of all motions, which are normal damp- 
ing, rolling, sliding and twisting. Figure |8] shows the dissipated 
energy for each mechanism. The solid lines represent the dissi- 
pated energies in the case without the normal damping and the 
dashed lines represent those in the case of k^ = 0.01. The dis- 
sipated energy in the case of k„ = 10 is indistinguishable from 
those in the case of k„ = 0.01, and thus we do not plot them. 
Note that the dissipation energy of the sliding force is less than 
lO"** erg, and thus it is not depicted in this figure. The dissipation 
by the rolling and twisting is almost the same in the cases with 
and without the normal damping. Thus, we confirm that the nor- 
mal damping does not affect the compression strength although it 




Fig. 8. Energy dissipation of each dissipation mechanism in [erg] 
against filling factor 0. The solid lines shows the result in the case 
without the normal damping and the dashed lines in the case of k„ = 
0.01 and The results in the case of k„ = 10 are not plotted because they 
are the same as those in the case of k„ = 0.01 and indistinguishable. 
The dissipation mechanisms are noimal damping, rolling, sliding and 
twisting. The dissipation energy by sliding motion is less than 10"' erg. 



dissipates the energy of the normal oscillations. Aside from the 
normal dissipation, the dominant dissipation mechanism is the 
rolling motion. This clearly shows that the static compression is 
determined by rolling motion of each connection, as mentioned 
in Section im Where > 0.1, the energy dissipation by twist- 
ing motion occurs. This is why Equation (l25T l is valid until the 
filling factor reaches 0. 1 as mentioned in Section lTTI In the high 
density region, where (f> > 0.1, another formulation is required 
but that is beyond the scope of this paper 

3.5. Dependence on the rolling energy 

We also investigate the dependence of the compression strength 
on ^crit- Since Ziion is proportional to ^cnt, we investigate the 
dependence on the rolling energy in this section. Figure |9]shows 
that the dependency on ^cnt- We vary ^cnt with 32, 16, 8, 4, 
and 2 A. The fixed parameters are = 16384, Cv = 3 x 10"^ 
and - 10"^. This result shows that the compression strength 
is almost the same in the low density region. This is because 
the periodic boundary creates the additional voids as discussed 
in Section im and thus we should not focus on the low density 
region. The lines in the case of ^cnt = 2,4, and 8 A are on their 
corresponding lines of Equation (l25t where (p < Q.l. The line in 
the case of ^cnt = 16 A has a little deviation and that in the case 
of ^crit = 32 A has a deviation from their corresponding lines of 
Equation dZST l. The reason why the lines in the case of ^cnt = 
16, 32 A deviate from the corresponding lines of Equation (l25t 
is that the dissipation energy is dominated not by rolling motion 
but by twisting motion as indicated in Figure [10] This figure 
shows that dissipated energy of each dissipation mechanism. We 
show the results of the cases with ^cnt = 8, 16, and 32 A. The 
normal damping is not contribute to the compression strength as 
discussed in Section 13.41 and thus we focus on the rolling and 
twisting motions. 

When ^crit < 8A, the dissipation energy is dominated by 
rolling motion. In the case of ^ciit - 32A, on the other hand, 
the dissipation energy is dominated by twisting motion. In the 
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case of ^crit - 16 A, the dissipation energy of rolling and twisting 
motion is comparable and thus this is the marginal case. Thus, 
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Fig. 9. Pressure P in [Pa] against filling factor with different critical 
displacement, fcTit- We put the same ten initial conditions varying ^crit 
with ^crit = 32, 16,8,4, and 2 A, respectively. Each line shows the 
average of ten runs. The other parameters are A' = 16384 Cv = 3 x 10"^, 
andfe„ = 10"^ 
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Fig. 10. Energy dissipation of each dissipation mechanism in [erg] 
against filling factor 0. Each panel represents the case of different <f;.,it, 
which are 32, 16, and 8 A, correspondent to Figure |9] 



the reason why Equation (l25l l is not valid when ^crit ^ 16A is 
that the twisting motion is the dominant mechanism to determine 
the compression strength. Therefore, we conclude that Equation 
is valid when fcrit ^ 8A. 



3.6. Fractal structure 

We also investigate how the fractal structure of the dust aggre- 
gate changes. Figure[TT]shows how many particles are inside the 
distance Hn for four snapshots. We select one run from the case 




Fig. 11. Number of particles inside the radius ri„ against normalized 
radius Hn/'o- This figure represents the fractal structure of the com- 
pressed aggregates in our simulation for various (p. We set a particle as 
an origin and count the number of particles inside r < Hn, where r is 
the distance from the origin to each particle's center. Then we count 
the same correlation of all particles as a n origin and take th eir average 
(similar figure of Figure.7 in the paper of lWada et al.l ( l2008i) ). Each line 
shows the result at the different time step. The solid thick lines repre- 
sents the structure of fractal dimension D = 2, and dashed lines repre- 
sent £) = 3 for each corresponding (p. The dotted line shows the number 
of particles in calculation. The region below this line corresponds to 
inside the periodic boundaries. 

with = 16384, Cy = 3 x 10"^ k„ = 10"-, and ^erit = 8 A. 
Each snapshot is when (p - 0.003,0.01,0.03 and 0.1, respec- 
tively. We take a particle as an origin and count the number of 
particles inside r < where r is the length from the origin. 
Then we set for all the other particles inside the computational 
region as an origin and take an average of them. We obtained 
the same trend in several runs in the cases of different shapes of 
initial aggregates. 

Note that we also count particles beyond the periodic bound- 
aries. In high Tin, oc r^^^ because copies over the periodic 
boundary distributed as fractal dimension of 3. Therefore, where 
N{r < Tin) > 16384, A^ must be A oc rK However, it is almost 
out of range of Figure [TT] The dotted line in Figure [TT] shows 
the number of paiticles in calculation, which is N - 16384. The 
results over this line is affected by the periodic boundary condi- 
tion and those below this line is in computational region. Thus, 
the results below the line represents the fractal structure inside 
the computational region and are not the artificial effect of the 
periodic boundary condition. 

Since the initial aggregate is a BCCA cluster, N is propor- 
tional to . In the case of = 0.003, which is equivalent to (f> of 
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the initial BCCA cluster, oc rr^ as shown in Figure [TT] When 
the fractal dimension is 3, can be written as 



N{r < nn) = 



Vo 



(29) 



where V(r < ri„) = (4/3)7rr^^. We also plot this equation as 
dashed lines for each in Figure [TT] Each dashed line has a 
good agreement in the large scale, while maintaining oc in 
small scale. 

Therefore, the structure evolution in the static compression 
is as follows. Initially, oc r^^ because the aggregate is a BCCA 
cluster As compression proceeds, the fractal dimension D be- 
comes 3 in a large scale while it is 2 in a small scale. The transit 
scale from D = 2 to D = 3 becomes smaller as compression pro- 
ceeds until D - 3 in any scale. This structure evolution means 
that the static compression reconstructs the aggregate first in a 
large scale with keeping the small scale BCCA structure. This 
is the reason why the rolling motion determines the compression 
strength, as discussed in Section |4] 



3.7. Silicate case : Comparison witli previous studies 

The compression strength has been investigated in the previous 
study (Seizingeret al. 2012). To investigate the connection of 
compression strength from the low density to the high density re- 
gion, we perf orm simulations i n the c ase of silicate with the same 
parameters of lSeizinger et al.l(l2012l) . Figure [T2lshows compres- 
sion in the case of silicate whose monomer size is 0.6 fim. The 
parameters are A^ = 16384, Cy = 3 x 10"^ and k„ = 0.01. The 
solid lines in Figure [T2 a) show the results of ten runs with dif- 
ferent initial aggregates and the thick solid line in Figure \T7\ h) 
shows their average. Using the rolling energy of silicate, which 
is £1011 = 1.42 X 10 erg, we also plot the line of Equation ( |25] | 
in Figure [T2l b). Since fo is given by 1.71 x 10"^ sec in the case 
of silicate aggregates, vj becomes 4.01 cm/s for <p = 10"^ with 
Cy-3 X 10"'. This Vd is larger than Cs etf (= 0.77 cm/s when 
4> = 0.01) for silicate aggregates, allowing the numerical results 
shown in Figure [12] to deviate from the line of Equation ((25) in 
the low region. When Vd = Cs.eff, (f> - 3.4 x 10"^, and there- 
fore, the compression strength should obey Equation ((25) when 
4> > 3.4x10^2 jnthe case of silicate, computational time is huge 
compared with ice particle cases. We take relatively high value 
of the boundary speed to save the computational time. Therefore, 
the result is deviate from Equation ((25) in the low density region 
because of the high velocity. In other words, the compression is 
not static in the low density region. In the high density region, 
on the other hand, the result is in good agreement with Equation 
((25), suggesting that Equation ((25) is applicable to aggregates 
consisting of silicate particles with different ro. 

To directly compare with previous studies. Figured?] shows 
the filling factor in linear scale agai nst pressure in log scal e. This 
figure corresponds to Figure 4 in (Seizinger et al.l ((2012[ ). The 
solid lines are our simulation results and the dashed Une is Equa- 
tion ( (25) in the low density region (0 < 0.1). T he dotted and 
solid lines are the result of Se izinger e ^^ ( 2012b and the fitting 
formula to experiments (iGiittler et al.ll2009 ). respectively. They 
performed similar A^-body simulations to ours but using a BPCA 
aggregate composed of silicate particles as an initial condition. 
The compression strength of our simulations has a good agree- 
ment with the same interaction model in Seizinger et al. (2012) 
with a little discrepancy; = 0.24 at P = 300 Pa in our si mula- 
tions and = 0.21 at P = 300 Pa in [Seizinger et al.l((2012[) . The 
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Fig. 13. The filling factor (f> against pressure P in [Pa]. This figure is 
same as Figure [T2] but plotted with linear scale of an d reversal of xy 
axis to compare with previous studies (see Figure 4 in ISeizinger et all 
The dotted line is the result of numerical simulations in the 
high density region ((f> > 0.1) in ISeizinger et al.l (120121) and th e thin 
solid line is the fitting formula proposed by iGuttler et al.l (l2009h . Our 
results consistently connect to the previous simulations in the high den- 
sity region. 



discrepancy, I3%m<p may be caused by the difference in the ini- 
tial aggreg ate or the pre s sure m easurement method. The fitting 
formula of (Guttler etaD(l2009l) suggests = 0.17 at P = 300 in 
the experiments. The discrepancy from our simulations is 29 % 
in (f>. In applicable uses of the static compression formula, we 
focus on obtaining ^ with a given P. 



4. Understanding the compression strength formula 

In this section, we analytically derive the compression strength 
and confirm Equation ((25) . First, we consider the structure of a 
fluffy aggregate in static compression in our simulations. As de- 
scribed in Section (273] we adopt the periodic boundary condition 
and put a BCCA cluster as the initial condition. This corresponds 
to a large aggregate which filled up with BCCA clusters three di- 
mensionally. As compression proceeds, the initial BCCA cluster 
is compressed but the aggregate keeps smaller BCCA structure 
as confirmed in Section 13.61 Therefore, the aggregate in static 
compression always consists of BCCA clusters in some scale 
and filled up with them. Figure [T4[ illustrates the aggregate in 
static compression. The enclosed lines depict BCCA clusters in 
a small scale. 

Next, we consider why the compression strength can be de- 
termined by the rolling energy. The internal mass density and 
the volume filling factor of the aggregate are equal to those of 
the BCCA clusters. Compression of the whole aggregate pro- 
ceeds by compression of each cluster Therefore, the compres- 
sion strength of the whole aggregate would be determined by 
BCCA clusters. The right panel of Figure [14] illustrates com- 
pression of one of BCCA clusters. The pressure on the BCCA 
cluster is exerted by neighbor clusters, which causes the com- 
pression of the BCCA cluster The BCCA cluster can be further 
divided into two smaller subclusters because BCCA clusters are 
created by cluster-cluster aggregation. A large void exists be- 
tween the two smaller clusters and they are connected with one 
connection of monomers in contact, represented by dashed line 
in the right panel of Figure [14] The compression of the BCCA 
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Fig. 12. Pressure P in [Pa] against filling factor (f>. This figure is same as Figure|4]but for the case of silicate particles (ro = 0.6/jm). 
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Fig. 14. Schematic drawing of compression of a dust aggregate consisting of a number of BCCA clusters. The left figure shows a dust 
aggregate consisting of many BCCA clusters and the BCCA clusters are distributed three dimensionally. Each enclosed line represents each 
region dominated by the BCCA clusters. The central figure is a BCCA cluster, receiving pressure from the next clusters. The BCCA cluster has a 
large void depicted in the central figure, and thus the void would be compressed, as expressed in the right figure. The required energy to compress 
the void is the energy to rotate the connection of monomers in contact. Therefore, the compression can be determined by the rolling motion of 
monomer connection on the connecting point of the subclusters. 



cluster occurs by crashing the large void, which requires only 
rolling of the monomers at the connection. 

Now, let us estimate the compression strength. In static com- 
pression, the aggregate is compressed by external pressure. Each 
BCCA cluster feels a similar pressure, P. Using the pressure, the 
force on the BCCA cluster is approximately given by 



Since the crashing of the large void is accompanied by rolling of 
a pair of monomers in contact, the work required for the crash - 
ing is given by so-c a lled th e rolling energy of monomers, £,oii 
dPominik & TielensI (Il997h or see Equation ([T) for its defini- 
tion). Therefore, the required force to compress the aggregate 



satisfies. 



P ■ /"BCCA ~ -Ell 



■oil- 



(31) 



Substituting Equation dSOl l, we further obtain the required pres- 
sure to compress the aggregate as 



(30) p 



PmW 



'BCCA 



(32) 



The radius of the BCCA clusters can be written by using the 
physical values of the whole aggregate. The internal density of 
the BCCA cluster is dependent on its radius. The BCCA cluster 
has the fractal dimension of 2, and its radius is approximately 
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given by tbcca = N^^^fo, where is the number of constituent 
monomers in the BCCA subcluster The internal density of the 
BCCA cluster is evaluated as 



NniQ I rficcA \ 



(33) 



'BCCA 



Using equations ( l32l i and ( l33T l, we finally obtain the required 
pressure (or the compression strength) as 



£1-011 / P 
rl \Po 



(34) 



This is the same as Equation dZSl ) obtained from our numerical 
simulations. 



5. Summary 

We investigated the static compression strength of highly porous 
dust aggregates, whose filling factor is lower than 0. 1 . We per- 
formed numerical A^-body simulations of static compression of 
highly porous dust aggregates. The initial dust aggregate is as- 
sumed to be a BCCA cluster The particle-particle interaction 
model is based on Dominik & Tielens (1997) and Wada et al. 
(l2007l) . We introduced a new method for compression. We 
adopted the periodic boundary condition in order to compress 
the dust aggregate uniformly and naturally. Because of the pe- 
riodic boundary condition, the dust aggregate in computational 
region represents a part of a large aggregate, and thus we could 
investigate the compression of a large aggregate. The periodic 
boundaries move toward the center and the distance between the 
boundaries becomes small. To measure the pressure of the ag- 
gregate, we adopted a similar manner used in molecular dynam- 
ics simulations. As a result of the numerical simulations, our 
main findings are as follows. 

- The relation between the compression pressure P and the fill- 
ing factor (f) can be written as 



P = 



^roll 



(35) 



where £1011 is the rolling energy of monomer particles and 
ro is the monomer radius. We defined the filling factor as 
(p - pi Pa, where p is the mass density of the whole aggre- 
gate, and po is the material mass density. Equation (l35T l is in- 
dependent of the numerical parameters; the number of parti- 
cles, the size of the initial BCCA cluster, the boundary speed, 
the normal damping force. We confirmed that Equation (l35t 
is applicable in different E^oW and r^. We also analytically 
confirmed Equation (l35T l. 

Equation ( |35] ) is valid where < 0.1 in the high density re- 
gion. In the low density region, we confirmed that Equation 
( |35] ) is valid for 4> > 10"^ in the case of = 16384. From 
the results of different initial sizes of the aggregates. Equa- 
tion ( l35T l is valid in the lower density region in the case of 
the larger aggregates. 

The initial BCCA cluster has a fractal dimension of 2 in 
the radius of the cluster, although the whole aggregate has 
a fractal dimension of 3 because of the periodic boundary. 
As compression proceeds, the fractal dimension inside the 
radius of the initial BCCA cluster becomes 3, while the frac- 
tal dimension in smaller scale keeps being 2. This means 
that the initial set up, which is that fractal dimension in large 



scale is 3 and that in small scale is 2, well reproduce the 
structure of a dust aggregate in static compression as a con- 
sequence. This also supports the fact that the compression 
strength is determined by BCCA structure in a small scale. 
- The static compression in the high density region (0 > 0.1) 
has been investigated in silicate case in previous studies 
( ISeizinger et al.ll2012l) . We performed the numerical sim- 
ulations in silicate case and confirmed that our results are 
consistent with that of previous studies in the high density 
region. 

The compression strength formula allows us to study how 
static compression affects the porosity evolution of dust aggre- 
gates in protoplanetary disks. In application to dust compression 
in protoplanetary disks, we use the compression strength for- 
mula to obtain with a given P. Moreover, the obtained com- 
pression strength would be applicable to SPH simulations of dust 
collisions. Such application of the static compression process is 
important future work. In this work, we did not study shear or 
tensile strengths. Those strengths are worth investigated in fu- 
ture work. 
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